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Abstract:  This  report  is  motivated  by  the  need  to  understand  the  possible  applications  of  two- 
dimensional  (2-D)  modeling  of  stochastic  processes  and  related  techniques  in  signal  processing  to 
problems  of  inversion  of  geophysical  data.  Here  we  first  provide  a  brief  overview  of  2-D  modeling 
and  signal  processing  techniques.  We  then  address  the  problem  of  passing  from  free  air  gravity 
anomaly  data  (FAG)  to  bathymetry.  We  present  an  approach  to  this  fundamental  inversion  prob¬ 
lem,  based  on  a  systematic  spatial  segmentation  of  the  2-D  data  followed  by  statistical  modeling 
within  the  segments.  We  also  suggest  possible  improvements  on  current  approaches  based  on 


transfer  function  methods. 
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1.  Introduction: 


It  Is  now  possible  to  obtain  very  reliable  mappings  of  the  earth’s  gravity  field 
using  air-borne  gravimeters.  Measurements  of  sea-level  using  radar  or  laser 
altimeters  on  satellites  also  provide  a  wealth  of  Information  concerning  the  geold. 
Typically  the  data  is  collected  along  tracks  followed  by  the  Instrument-bearing 
platform.  The  tracks  provide  a  grid  like  pattern  of  data  to  be  used  for  Inversion 
purposes,  e.g.  Inverting  FAG  to  bathymetry.  Current  practice  In  this  field  Is 
based  on  modeling  and  signal  processing  methods  that  are  Intrinsically  one 
dimensional  (along  the  data-collectlon  tracks)  In  nature.  To  make  effective  use  of 
such  techniques,  It  Is  also  customary  to  assume  certain  directional  uniformity  pro¬ 
perties  perpendicular  to  the  tracks.  From  the  point  of  view  geophysics,  such 
assumptions  may  be  unnatural.  Also,  processing  along  tracks  often  Implies  only 
partial  use  of  the  Information  at  hand.  These  and  other  reasons,  motivated  us  to 
Investigate  the  possibility  of  using  the  methods  of  two  dimensional  stochastic 
modeling  for  the  purposes  geophysical  Inversion. 

In  the  present  report,  we  begin  by  giving  a  brief  overview  of  certain  aspects 
of  the  theory  of  2-D  random  fields  In  section  2.  We  focus  primarily  on  certain 
recent  developments  In  2-D  autoregressive  moving  average  (ARMA)  type  models. 
An  Important  Idea  In  our  work  Is  that  when  such  models  are  constructed  for 
FAG  data,  they  will  be  valid  only  over  limited  regions  of  the  earth,  In  part  deter¬ 
mined  by  the  ’largest  wavelength’  appropriate  for  the  data.  It  Is  therefore  neces¬ 
sary  to  properly  segment  the  FAG  data  Into  various  spatial  regions  and  construct 
models  within  the  regions.  Our  approach  to  segmentation  Is  based  on  detecting 
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edges  across  which  large  changes  In  FAG  data  arise.  We  discuss  a  class  of  tech¬ 
niques  for  detecting  edges  In  FAG  data  In  section  4.  (The  methods  may  be  of 
Independent  geophysical  Interest  than  for  Inversion.) 

In  section  4.1,  we  present  our  approach  to  the  use  of  two  dimensional  model¬ 
ing  techniques  for  the  problem  of  Inverting  from  FAG  data  to  bathymetry.  We 
give  details  of  our  estimation  technique  and  segmentation  technique  In  sections 
4.2  and  4.3.  In  section  4.4,  we  present  a  brief  discussion  of  some  recent  advances 
In  the  theory  of  deconvolution.  We  think  that  these  new  results  could  be  of  con¬ 
siderable  use  In  modeling  and  Inversion  based  on  transfer  functions. 

The  final  section  5  contains  proposals  for  future  work.  Among  other  things, 
we  believe  that  detailed  studies  of  sea  bottom  roughness  models  may  help  make 


our  statistical  techniques  more  reliable. 
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2.  Parametric  Models  for  2-D  Random  Fields: 

A  basic  problem  associated  to  an  array  of  (random)  data  on  a  two- 
dimensional  “grid”  Is  that  of  reconstructing  the  underlying  signal  corresponding 
to  the  data.  Based  on  their  success  In  treating  one  dimensional  signals, 
parametric  models  are  highly  desirable  In  this  Instance  also.  Two  dimensional 
analogues  of  ARMA  models  can  be  Introduced  for  this  purpose.  We  do  so  follow¬ 
ing  a  preliminary  discussion  of  random  fields. 

2.1.  Random  Fields 

A  stochastic  process  over  T  Is  a  family  of  random  variables  Indexed  by  a 
parameter  set  T.  We  denote  this  family  as  x  =  {Xt :  t  e  T}.  Here  we  assume 
that  the  random  variables  Xt  take  vallues  In  a  vector  space  -  typically  <XN .  If  N 
=  1,  It  Is  customary  to  refer  to  x  as  a  univariate  stochastic  process. 

Usually,  the  Index  set  T  denotes  time,  l.e.  T  /=  ]R  or  Z .  However  If  T  Is  a 
multidimensional  space,  we  say  that  x  Is  a  random  field.  For  our  purposes,  we 
will  mostly  restrict  attention  to  T  =  nt2  (random  field  on  a  plane)  or  T  =  Z2 
(random  field  on  a  planar  grid  or  lattice). 

Many  Important  concepts  from  the  theory  of  stochastic  processes  over  K  (or 
Z)  generalize  to  random  variables  over  Rn  (or  Zn). 

We  say  that  a  univariate  random  field  x  over  Rn  (or  Zn)  taking  values  In  CE 
Is  homogeneous  If, 

E(Xt)=n  -VteR”  (or  zn). 


and 
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E(Xt+TX, t  +r)  =  E  (Xt  Xf  ),  \/t ,  «  ,r£R‘  (or  Zn). 

We  denote  a s  R(t  -  tf  )  the  correlation  E  {(Xt  -  E(Xt))(Xe  -  E  (Xf  ))}  of  a  given 

homogeneous  random  field.  It  Is  Important  to  note  that  R(t-f  )  Is  non-negative 
definite,  l.e.,  given  a  set  {a„  .  .  .  ,  am }  of  complex  numbers,  m  an  arbitrary  posi¬ 
tive  Integer, 


E  E  a<aj  R  (((  -  </)  >° 

i=l  ;=  1 

for  any  choice  }  C  ir"  (or  Zn ). 


It  follows  by  a  famous  theorem  of  Bochner  and  Herglotz  that  we  have  a 
(spectral)  representation, 


R(t)=  J  ei<v,t>  dF (u). 

H”  (2.1.1) 

or 
j>  n 

(Here  Tn  denotes  the  n-torus).  F  Is  a  measure  known  as  the  spectral  measure. 
When  It  Is  possible  to  write, 


dF(v)  =  S(u)d  o, 

we  say  that  S  (ifi  Is  the  spectral  density  of  the  given  random  field. 

Homogeneous  random  fields  are  singled  out  for  modeling  data  sets  that 
appear  to  possess  a  translation  Invariance  property  In  their  correlation.  A  further 
special  class  of  random  fields  of  Importance  Is  the  class  of  homogeneous  and  iso¬ 
tropic  random  fields  satisfying; 

E  (Xt )  =  n  =  0  \/t  E  IR"  (or  Zn), 

and 


E(Xt,X?  )  =  E(XP{t)XP{v  ))  \/t,V  €  JRn  (or  Zn), 

and  the  map  t^-P(t)  denotes  a  rigid  motion  In  E"  (or  Zn ).  Recall  that  a  rigid 
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motion  In  R"  Is  a  combination  of  rotation  and  translation.  Thus,  for  a  planar, 
homogeneous  and  Isotropic  random  field  to  be  a  good  model  for  a  given  data  set, 
the  latter  should  have  both  translational  and  rotational  Invariance  properties  In 
Its  statistics. 

It  follows  that 


E(XtX(  )  =  i?(  |  |  t  -  H  |  |)  (2.1.2) 

Is  a  function  only  of  the  Euclidean  distance  |  |  t  -  V  |  j  (or  the  lattice  dis¬ 
tance).  A  very  nice  theorem  analogous  to  Bochner-Herglotz  holds. 

Theorem: 

A  function  R  (n),  0<i/<oo  Is  the  correlation  function  of  an  Isotropic  and 
homogeneous  random  field  over  IR”  Iff, 


r  J f»-2)/2(X0 

R  (r  )  =  /  ■  {—Rr-:.,dF 0(X). 
o  (Xr)<n-2)/2 


(For  n=2,  R  (r  )  =  f  J0(\r  )  dF0(\)) 


(2.1.3) 


Here  J0—  Bessel  function  of  order  0,  etc. 

Remark  2.1.1:  This  theorem  tells  us  how  to  pass  correctly  from  correlation  func¬ 
tions  to  spectral  density  functions.  The  ordinary  Fourier  transform  Is  not  the 
right  thing  to  use  for  homogeneous  and  Isotropic  random  fields. 

Remark  2.1.2:  For  data  In  2  dimensions,  a  plot  of  the  correlation  function  should 
reveal  radial  symmetry  If  present  and  suggest  what  class  of  random  field  models 
to  use.  We  will  In  any  case  always  make  the  assumption  that  homogeneity  Is 
valid  for  FAG  data  if  the  data  is  restricted  to  a  reasonably  sized  segment.  Pre¬ 
cisely  how  to  determine  the  segmentation  will  be  one  of  our  major  concerns. 
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Before  we  close  this  section,  we  note  that  the  property  of  Markovlanness 
which  plays  a  key  role  In  connection  with  stochastic  processes  over  ir  has  a  con¬ 
siderably  more  complex  theory  In  the  case  random  fields.  One  definition  based  on 
conditional  Independence  has  the  associated  picture: 


{Xz :  zeiR2}  Is  Markov  If  for  any  region  D  -  as  above,  X2_and  X2+  are  Independent, 
conditioned  on  {X2 :  z  e  dD  }.  We  shall  see  below  In  section  (2.3)  that  a  different 
notion  of  Markovlanness  Is  needed. 

Since  we  are  primarily  concerned  with  grldded  FAG  data,  we  shall  confine 
all  our  discussions  on  modeling  to  2-D  random  fields  over  the  planar  lattice. 

2.2.  Models  for  2-D  Discrete  Random  Fields: 

From  2-D  linear  system  theory,  we  have  the  notion  of  an  Input-output  sys¬ 
tem  as  a  discrete  convolution; 

Y (n  vn 2)  =  (/ -H  *  X){nltn2) 

00  00 

=  S  E  h(n1- m„  n2- m2)X(m1,  m8)  (2.2.1) 

m  j=  -oo  »i  2~  ~  00 

where  X(v)  Is  an  input  and  F(v)  Is  the  corresponding  output  and  h  (•,  )  denotes 
the  impulse  response  or  weighting  pattern. 
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Tlie  system  (2.2.1)  Is  causal  If, 

h  (n  un2)  =  0  Vw1(  n2<0.  (2.2.2) 

The  weighting  pattern  h  {■,■)  Is  said  to  be  separable  If 

h  (nl,n2)  =  /  (nx)g  (rt2)  \/nun2eZ.  (2.2.3) 

Much  of  2-D  realization  theory,  [Bose  1082,  Multidimensional  System  Theory]  Is 

concerned  with  conditions  on  weighting  patterns  satisfying  (2.2.2)  that  can  be 

realized  as,  difference  equations  of  the  form, 

u,  m2 

E  E  f(ni-m1,«rm,)  (2.2.4) 

TO  J— 0  TO  2=0 

lx  l2 

—  E  E  Aj,<2  1  -  ^ i’  n2  _ 

L=°  (2=° 

If  T(v)  and  X(v)  are  random  fields  and  If  we  let  X(v)  be  a  discrete  2-D  white 
noise  process,  then  (2.2.4)  becomes  the  analogue  of  the  notion  of  ARMA  model 
widely  used  In  time  series  analysis  [Box  Jenkins  1971  Time  Series]. 

We  say  that  H  Is  stable  If  In  (2.2.1)  we  get  bounded  outputs  for  bounded 
Inputs.  We  note, 

Theorem,  (Stability): 

00  00  -  n  -  n 

Let  H(z1,z2)=  Yj  S  h(n1,n2)z1  1 z2  2  and  assume  that  H  Is  casual 

n  j=  -  00  n  2=  -  00 

and  of  the  form,  H(zuz2)  =  ^ZuZP .  Then  H  Is  stable  If  It  Is  analytic  In 

B(zx,z2) 

I  Z,  |  >1,  |  Z*\  >1. 

While  the  subject  of  causal  2-D  models  Is  of  Independent  Interest  there  Is  no 
compelling  reason  for  us  to  work  with  such  models  In  the  present  context  of  geo¬ 
physical  fields.  One  would  however  expect  geophysical  data  to  satisfy  certain 
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relations  between  field  values  associated  to  points  In  a  given  lattice  neighborhood 
of  a  lattice  point.  Thus,  we  are  In  a  situation  where  noncausal  2-D  models  gen¬ 
eralizing  (2.2.4)  are  needed.  One  class  of  such  models  Is  given  by  the  following: 

Xr  ,t  =  X;  ai(Xr-i,t  +  Xr+lV)  + 

i  =1 

4  Yj  b,  (Xr,i-i  4  Xr  ,+J  )  +  WT  ,  (2.2.5) 

j  =i 

Here  {Wr , ;  -  co<r  <oo  ,-o o<s  <oo}  Is  the  2-D  analog  of  discrete  time  white  noise, 
satisfying 

E(Wri,)  =  o 

cov(Wrt,  ,Wy  y  )  =  S(r-rl  ,  s -s'  )q,  q  >0 
(8(i  ,j )  Is  the  usual  Kronecker  symbol).  Models  of  this  type  are  Investigated  In 

the  elegant  volume  [Bartlett  1976,  Spatial  Patterns].  In  the  special  case, 

-^r+l,  »+l  —  al-^r,i+l  4  (2.2.6) 

j 

-a  ia2Xr,,  +  Wr[, 

we  get  an  exponential  correlation  of  the  form, 

Pxx(*.l)  =  E  (Xr  Xr  +i,,+i ) 

=  <x2exp[  -  c !  I  i  I  -  c  2  I  y  I  ]  «!,  cs>0.  (2.2.7) 

The  parameter  cv  c2  and  a2  In  the  correlation  function  are  related  to  the  param¬ 
eters  of  the  casual  model  (2.2.6)  by  the  following  formulas: 

a !  =  exp(-c  ,) 

a2  =  exp(-c2)  (2.2.8) 

q  =  cr2(i  -  a2){  1  -  a2 ). 

From  (2.2.8),  It  Is  clear  that  one  can  estimate  the  parameters  of  the  model  (2.2.6), 
If  Indeed  the  data  justifies  a  correlation  function  of  the  form  (2.2.7). 
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Non-causal  models  of  the  form  (2.2.5),  (2.2.6)  and  their  variants  are  of 
Interest  to  us  as  candidates  for  geophysical  problems.  There  are  two  fields  of 
Interest  to  us:  (a)  the  free  air  gravity  (FAG)  field  and  (b)  the  bathymetry  field. 
In  what  follows,  we  shall  be  concerned  with  Imposing  the  proper  constraints  on 
noncausal  models  In  order  to  be  useful  for  the  solution  of  geophysical  problems. 

2.3.  Periodic  Random  Fields: 

The  random  fields  of  section  (2.2)  were  assumed  to  be  defined  over  the  entire 
planar  lattice.  This  conflicts  with  two  Important  factors: 

o  for  a  geophysical  field  such  as  FAG  data  or  bathymetry,  one  cannot  reason¬ 
ably  expect  spatial  homgenelty  properties  outside  a  certain  window  or  seg¬ 
ment;  different  parts  of  the  world  should  obey  different  local  regularities 
(models). 

o  practical  limitations  of  data  collection  confine  one  to  such  a  window  to  begin 
with. 

Thus  we  are  forced  to  Impose  additional  window  constraints.  It  Is  quite  obvious 
that  these  considerations  are  not  unique  to  2-D  models.  In  fact  the  above 
remarks  apply  equally  well  to  1-D  models.  In  his  book  on  time  series,  Hannan 
[Hannan  1960,  Time  Series]  considered  a  class  of  processes  that  are  periodic,  l.e. 
satslfy  a  constraint  of  the  form, 

n+«m=n  Vk.mez  (2.3.1) 

(N  =  window  size) 

In  the  present  context,  one  can  apply  the  same  Idea  to  a  noncausal  random  field 
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model,  namely  to  define  the  field  simply  by  periodic  extension  outside  a  window 
of  size  N1xN2.  As  we  shall  see  below,  this  has  Important  consequences  for  com¬ 
putations  as  well.  In  several  recent  papers,  Kashyap  and  Chellappa  have 
analyzed  an  Important  special  class  of  noncausal  2-D  periodic  random  field 
models  [Kashyap  1981,  Image  Modeling]  [Chellappa,  Kashyap  1979  Decision  Rules] 
[Chellappa  1980  Spatial  Autoregressions].  For  a  variety  of  reasons  that  will  be 
apparent  In  section  4,  we  adopt  precisely  this  class  of  models  as  candidate  models 
for  geophysical  fields. 

Let  a  2-D  discrete  real  random  field  { Yitj :  i  ,je  Z }  be  given  by  an  autoregres¬ 
sive  model  of  the  form, 


Y*.i  +  Yi  +  qkl,  i  +  qk2  (2.3.1) 

here  Q  :={(qk\  qkz):  k  =i,2,...,m  }  Is  a  subset  of  Z  xZ  Identifying  the  appropriate 
neighborhood  of  lattice  points  associated  to  a  given  one.  The  9{  are  scalar 
weights,  Uii}-  Is  a  discrete  2-D  random  field  corresponding  to  Gaussian  white 
noise;  E  [Ify  ]  =  o  and  Var  [U{j  ]  =  l. 

The  model  (2.3.1)  Is  taken  together  with  a  periodicity  condition, 


Yi  +  kN  j.y  +  =  Yi.j  Vi  ,j  ,k  ,1  €  Z.  (2.3.2) 

Here  Nk  and  N2  are  positive  Integers  determining  the  window  or  segment  of 

homogeneous  behavior.  We  shall  be  concerned  with  modeling  using  only  data 

from  within  the  window 


W:={(*\j):  i<t  <Nk;  1  <j  <7V2}. 


(2.3.3) 
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A  typical  example  of  Q  might  be 

Q  :={( — 1,0),  (0,-1),  (1,0),  (0,1),  (1,1)}.  (2.3.4) 

The  corresponding  neighborhood  can  represented  as  In  Fig.  2.3.1. 

(i  +  1  ,  j  +•  ^  ) 


(  iL- 1  .  i  }  :  •  -  ■  -•**  C 1  /  j  J 

C *), 


(jt ,  j  -  4- ) 

Figure  2.3.1 

In  section  4  below  we  analyze  the  class  of  models  (2.3.1)  -  (2.3.2)  and  give 
formulas  for  estimation  of  parameters.  It  should  be  borne  In  mind  that  the 
parameters  of  the  model  consist  of  the  tuple  (Q  ,  6,  Nlr  N 2,  0)  where  6  =  (Ou...,0m  )T . 

In  general  the  models  (2.3.1)  -  (2.3.2)  do  not  satisfy  the  usual  Markov  pro¬ 
perty,  l.e.  Prob  density  {Yitj  |  YliP;  (/,p)=^(* ,/)} 


T^Prob  density  {  Y,}  \  Ylf); 


(2.3.5) 
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l  =  i  +  qkl 
P  =  j  +  Qk2 
(9*1.  ?*2)£Q  } 

However  (2.3.5)  Is  Indeed  an  equality  If  Q  In  (2.3.5)  Is  replaced  by  some  Q  XZ>Q  . 

3.  Current  Practice  in  Geophysical  Inversion: 

The  main  problem  of  gravity  Interpretation  Is  that  of  Inferring  a  plausible 
subsurface  density  distribution  (or  subsurface  body)  from  surface  observations  of 
the  free  air  gravity  anomaly  (FAG).  More  specifically,  based  on  (FAG)  data  over 
the  ocean,  we  seek  to  Infer  the  bathymetry.  There  are  a  variety  of  sources  of 
nonuniqueness  In  the  inversion  process.  One  example  Is  that  of  density  distribu¬ 
tions  of  rapidly  oscillating  positive  and  negative  densities  at  depth.  Also  more 
seriously,  geophysical  processes  such  as  isostatic  compensation  lead  to  nonunique¬ 
ness.  We  make  some  remarks  about  current  practice  In  this  area. 

>  / 

3.1.  Discrete  Inverse  Theory: 

At  present,  two  basic  approaches  are  used: 

(a)  model  the  gravity  anomaly  by  means  of  of  one  or  more  constant- 
density  bodies  with  variable  geometry,  such  as  a  sphere,  cylinder,  prism, 
polygon,  a  set  of  prisms  to  model  an  Interface  etc.  The  geometry  of  the 
disturbing  body  Is  Iteratively  adjusted  from  a  starting  point; 

(b)  fix  the  geometry  of  the  model  (for  example  a  regular  array  of  Identical 
rectangular  or  square  blocks)  and  allow  the  densities  of  the  blocks  to  vary. 

In  connection  with  approach  (a),  both  spatial  and  frequency  domain  methods 


Figure  3.1  (Last  a  Kublk  1983) 

d  and  h  are  the  horizontal  and  vertical  dimensions  of  an  elementry  block.  Den¬ 
sity  of  j  th  block  is  v j . 
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have  been  used  [Oldenburg  1974  gravity  anomalies],  [Bhattacharya  Leu  1975 
gravity  and  magnetic  anomalies],  [Pedersen  1978  statistical  potential  fields], 
[Pedersen  1979  Wavenumber  domain]. 

In  connection  with  approach  (b)  there  are  substantial  difficulties  In  reducing 
the  nonuniqueness.  However,  recent  efforts  In  this  direction  Include  the  use  of 
quadratic  programming  and  monotonlclty  assumptions  on  the  density  [Fisher 
Howard  1980  Gravity  Interpretation  quadratic  programming]  and 
“compactlflcatlon,”  or  minimizing  the  volume  of  the  causative  body  [Last  Kublk 
1983  Compact  gravity].  We  recall  here  the  basic  model  of  [Last  Kublk  1983 
Compact  gravity].- 

Consider  Figure  3.1.  The  data  Is  collected  along  a  track.  The  geometry  Is 
fixed,  but  the  density  Is  allowed  to  vary  from  block  to  block.  The  model  for 
gravity  at  i  th  data  point  Is 

M 

9i  =  E  aijvi  +  e<  *  =  1,2,.. ..IV  (3.1.1) 

i  =i 

where  u;  —  density  of  j  th  block,  e,-  =  noise  associated  to  » th  measurement  and, 
ajj  —  matrix  element  representing  the  Influence  If  the  j  th  block  on  the  *th  grav¬ 
ity  value.  From  standard  arguments  [Garland  1965  Earth’s  Shape  and  Gravity] 
we  know, 

a, 7  =  27 [(*,-  -  Xj  +  d  /2)log(— —) 

r1ri 

+  d  log(— )  (3.1.2) 


-  {zj  +  h/2)(64  -  92) 
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+  (2Tf  -  A/2)(tfs  -  0,)] 

where, 


2 

ri 

=  (*y 

- 

h/2f 

+ 

(*f  -  xj 

+ 

d/2f. 

r  2 

=  (*/ 

+ 

h/2 ? 

+ 

{.X i  ~  X3 

+ 

d/2?, 

r  32 

=  (Z) 

- 

h/2? 

+ 

{Xi  -  Xj 

- 

d  /2?, 

r<2 

=  (zj 

+ 

h/2? 

+ 

(x{  -  Xj 

- 

d  /2 ?, 

—  tan" 

-  Xj 

+ 

d  /2)/{Zj 

- 

h/2), 

02  = 

=  tan- 

-  xi 

+ 

d/2)/(zj 

+ 

h/2), 

^3 

=  tan' 

~\*i 

-  xi 

- 

d/2)/(zj 

- 

h/2), 

0 4 

=  tan 

'  -  xi 

- 

d  /2)(Zj 

+ 

h/2). 

Here  7  Is  the  universal  gravitational  constant. 

It  follows  that  one  can  set  up  the  problem  as  that  of  solving  for  a  vector  V 
of  densities  satlsyflng, 

G  =  AV  +  E  ,  1  (3.1.4) 

where  G  Is  the  vector  of  gravity  data.  One  typically  solves  (3.1.4)  using  least 

squares.  There  are  various  approaches  to  regularizing  the  problem  Including  the 

compactlflcatlon  Idea  of  Last  and  Kublk.  It  Is  customary  to  use  weighted  least 

squares. 

More  recently  certain  Implementations  of  one  dimensional  track  Inversion 
have  become  available  that  are  Iterative  but  not  based  on  least  squares.  We  dis¬ 


cuss  one  below. 
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3.2.  The  Algorithm  of  Lipman: 

This  particular  algorithm  has  been  Implemented  and  tested  and  Is  discussed 
In  a  report  [Rosenblum  1983  gravity  to  bathymetry  Inversion].  Since  the  details 
are  available  In  the  report,  we  limit  ourselves  to  highlighting  the  main  features  of 
the  algorithm: 

(a)  the  method  Is  Iterative  starting  from  an  Initial  mean  bathymetry  such 
that  the  corresponding  FAG  ==  o; 

(b)  It  adjusts  the  bathymetry  by  adding  or  subtracting  mass  (via  prisms) 
at  each  point; 

(c)  It  assumes  uniform  density  along  a  track; 

(d)  the  Iteration  G  based  on  attempting  to  bring  the  computed  FAG  at 
each  point  to  some  value  not  further  than  a  tolerance  e  (In  mgals)  from  the 
corresponding  FAG; 

! 

(e)  It  Includes  a  constraint  on  slope  to  avoid  spiking; 

(f)  the  realization  of  step  (d)  above  Is  based  on  two  elements  -  a  decision 
rule  that  determines  whether  a  region  which  violates  the  predefined  toler¬ 
ance  6  In  (d)  above  Is  ‘fixable’  -  l.e.  Improved  to  within  tolerance  without 
adversely  affecting  points  outside  this  region  but  ‘nearby’. 

(g)  once  a  region  has  been  found  to  be  fixable,  new  prismatic  masses  are 
added  using  a  routine  called  SLABADD. 

The  above  algorithm  does  not  explicitly  use  a  goodness-of-flt  measure  as  In 
least  squares.  In  Its  current  Implementation,  the  algorithm  requires  the  following 


user-defined  Input  parameters: 
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(a)  the  FAG  data  point  spacing  In  the  data  file  (In  km); 

(b)  the  average  depth  along  the  track  (In  km); 

(c)  the  average  density  (assumed  constant); 

(d)  the  maximum  error  allowed  In  fitting  the  estimated  FAG  to  the  meas¬ 
ured  FAG  (In  mgals,  usually  2-4  mgals); 

(e)  the  bathymetric  stepslze  allowed  (In  kms,  about  1/5  the  data  point 
spacing  In  (a)  above); 

A  certain  amount  of  user-sophistication  appears  to  be  necessary  In  order  to 
ensure  adequate  performance  of  the  algorithm.  In  the  next  sub-section,  we  make 
some  comments  on  current  practice  In  order  to  set  the  context  for  our  proposed 
approach. 

3.3.  Some  Comments  on  Current  Inversion  Techniques: 

To  the  best  of  our  knowledge,  all  currently  used  techniques  are  based  on 
processing  data  along  tracks,  l.e.  these  are  1-D  methods.  This  Implies  that  for 
the  least  squares  methods  of  section  3.1.,  one  has  to  use  uniform  prisms  that  are 
of  Infinite  extent  perpendicular  to  the  tracks.  In  the  regions  of  the  earth  where 
such  directional  homogeneity  Is  unavailable,  the  use  of  1-D  methods  Is  at  best  a 
rough  and  ready  approximation.  However,  direct  extensions  of  the  methods  of 
section  (3.1)  to  2-D  data  lead  to  considerable  computational  problems.  Formulas 
analogous  to  (3.1.2)  and  (3.1.3)  can  be  written  down,  but  are  difficult  to  analyze 
for  the  purposes  of  determining  whether  or  not  “Influence  matrices”  such  as  A  In 
(3.1.4)  are  approximated  by  sparse  matrices.  Such  considerations  are  essential  In 
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order  to  make  the  problem  computationally  tractable. 

We  propose  that  a  different  approach  based  on  discrete  2-D  random  field 
models  be  considered  for  the  purposes  of  Inversion.  The  details  are  given  In  the 
next  section.  The  main  point  of  our  approach  Is  that  we  propose  model  building 
based  on  sound  statistical  principles.  The  class  of  models  that  we  consider  Is  that 
of  periodic  random  fields  Introduced  In  (2.3). 

4.  Geophysical  Modeling  and  Inversion  Using  2-D  Discrete  Random  Fields: 

In  this  section,  we  present  our  algorithm  for  modeling  and  Inversion  of  2-D 
random  fields  associated  to  the  problem  of  Inversion  from  FAG  to  bathymetry. 
The  basic  random  fields  for  us  are  the  univariate  fields  of  free  air  gravity  and 
bathymetry,  denoted  as  Gi}-  and  Bi}-  respectively.  For  computational  reasons,  we 
restrict  ourselves  to  univariate  techniques.  The  basic  steps  are  given  In  section 
4.1.  The  details  of  the  estimation  procedure  (following  Kashyap)  are  given  In  4.2. 
The  procedure  In  (4.1)  assumes  that  a  prior  segmentation  of  the  field  In  question 
has  been  given.  In  order  to  obtain  such  a  segmentation  or  windowing  we  discuss 
In  section  (4.3)  a  method  that  has  Its  origins  In  Image  processing.  This  Is  a  deter¬ 
ministic  two-stage  technique  Involving  first  a  convolution  with  the  gravity  field 
by  a  suitable  kernel  followed  by  nonlinear  processing.  In  comparison  with  a  sta¬ 
tistical  technique  proposed  by  Kashyap,  this  method  appears  to  be  much  less 
computation-intensive  and  has  performed  well  In  experiments  on  real  Images. 
Our  Idea  Is  to  apply  this  technique  to  the  gravity  anomaly  data  and  segment  It 


before  doing  any  detailed  modeling. 
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The  last  section  of  this  chapter  Is  denoted  to  some  recent  advances  In  decon¬ 
volution  techniques.  We  believe  that  these  new  tools  could  be  of  great  utility  In 
Inversion  algorithms  based  on  transfer  function  models. 

4.1.  The  Modeling  and  Inversion  Algorithm: 

We  use  the  notation  of  section  (2.3).  Recall,  from  that  section  that  a 
periodic  random  field  model  obeys  a  law  of  the  form: 

1 

m  — 

Yu  +  £  »k  Yt  i  ,  ,  =  /?2  U{j  1  <i  <NU  i <j  <N2, 

k=  1  *  '  * 

Yi+kNvi+tNt=  Yy  V«  ,j  ,k  ,1  ,E  Z.  (4.1.1) 

Here  N1  and  N  2  determine  the  window  size  and 

Q  =  {(qk\  qk2)  k  —  }c  Z  xZ  determines  neighborhood  dependence.  The 

field  { U(j }  Is  a  standard  discrete  Gaussian  white  noise  field.  We  find  It  con¬ 
venient  to  adapt  the  well-known  polynomial  model  notation  as  follows: 

Let  and  «2  denote  the  following  shift  operators, 

5 1  Yj  fj  =  L  _i,  j , 

s2Yij  =  Y{',+1  .  (4.1.2) 

Then,  the  difference  equations  in  (4.1.1)  take  the  form, 

Yi j  +  £  ekSl  k  s/*  l<i<Nul<j<N2  (4.1.3) 

k  =1 

or  more  compactly, 

©(*>.  -aW,  =^UiS,  C4'1’4) 


where 
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m  _■  1  2 

©(*i.  « a)  =  l  +  E  ^  s  i  *  «2  *  •  (4.1.5) 

A?  =i 

Each  such  model  Is  parametrized  by  the  t.uple,  (Q,  6,  Nu  Na,  0).  We  are  now 
ready  to  state  our  algorithms. 

Modeling  Algorithm:  (Known  G,y  and  Bi}- ). 

Step  1:  Choose  a  window/segment  of  size  NixN2  and  construct  a  bathymetry 
model 


mB  1  2 

(l  +  E  )Bii  ='f¥v?j  (4.1.6) 

k  =1 

Step  2:  Using  the  same  window/segment  construct  a  gravity  model 

mG  1  _  2 

(1  +  E  )G<i  =  WV#  (4.1.7) 

k  =1 

Step  3:  Investigate  whether  the  residual  field , 

Wfj  =  -  \fWu§  +  VjPUif 

passes  a  whiteness  test,  e.g.  by  explicitly  estimating  the  correlation  function  of 
W{j .  If  It  does  not  then  construct  a  residual  model  of  the  form, 

mW  _  1  _  2 

(1  +  E  W*  s2k  )Wii  =  sfWuis  (4.1.8) 

k  =1 

where  U{]-  Is  a  standard  white  noise  held. 

Step  4-'  We  now  obtain  a  model  relating  gravity  and  bathymetry, 


ma  1  _0  2 

(1  +  E  *  «9  *  )<?*■ 

it  =1 

1  -,2 

=  (1  +  E  «a  *  )%  + 

A?  =i 

mVK  1  2 

+  (1  +  E  *2~?‘  r1'/^TrC1-y  . 


(4.1.9) 
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Clearly,  this  can  be  written  In  the  form, 


©wr(Si.S2)‘©G(*i.s8)G,fy  =  ew(s1,s2)eB(s1,s2)Bi]- 

+  \fWUij  (4.1.10) 

The  model  Is  parametrized  by  the  following  neighborhood  sets  QB ,  Qa ,  Qw ,  the 
weight  vectors  9B ,  9a  ,  9W  and  the  window  size  N ^XN2. 

We  can  write  (4.1.10)  compactly  In  the  form, 

M (a  1(*  2)  Gij  =  N (s  „s  2)B,y  +  \Zj3wUij,  (4.1.11) 

where, 


N (5  1^2)  —  ^  (5  1 2)*  ®  (s  i>®  2) 

m(u) 

E 

k  —l 


1  +  E 


and 


M  (s  vs  a)  =  ©  w  {s  us  2)©°  ( s  us  2) 


(4.1.12) 


m  M 

=  1  +  E  Vk  » 

k  =1 


t* 


(4.1.13) 


where  m  (p)  and  m  {v)  are  positive  Integers  determined  by  mw ,  mB  and  m 


Suppose  we  are  now  presented  with  gravity  data  from  a  region  of  the  earth 
which  Is  geophysically  similar  to  the  region  from  which  the  data  for  gravity  and 
bathymetry  used  In  constructing  the  model  (4.1.11)  was  obtained.  It  Is  now 
required  to  Infer  the  bathymetry  for  this  new  region.  We  simply  postulate  that 

IV 

the  model  (4.1.11)  applies  In  the  new  region  as  well.  If1' this  new  region,  the  Inver¬ 
sion  process  has  thus  two  steps. 

is 

Inversion  from  Gi}-  /B,y  : 
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Step  1:  Postulate  a  model  of  the  form, 

M(Slls2)G{j  =  +  VjF'Uij 

with  known  polynomials  M  and  N  and  coefficient  pw . 

Step  2:  Given  gravity  data  G{j  In  the  new  region  solve  for  Bi}- ,  (e.g.  by  applying 
least  squares  to  the  above  model  or  by  constrained  least  squares  If  some  partial 
bathymetry  constraints  are  at  hand). 

From  the  discussion  above,  It  should  be  apparent  that  the  key  modeling  pro¬ 
cess  Involves  constructing  three  periodic  random  field  models  associated  to  the 
poynomlals,  6s(s,s2)  and  e"r(«1)os),  all  of  the  same  generic  form 

(4.1.1)-(4.1.2).  In  the  next  section  we  explain  the  basic  steps  In  this  estimation 
problem. 


4.2.  Estimation  of  a  Generic  Univariate  Autoregressive  Periodic  Random  Fields 

In  this  section  we  will  continue  to  use  the  notation  of  section  (2.3).  First  we 
recall  the  notion  of  a  clrculant  matrix,  l.e.  a  matrix  of  the  form 


circulant  )  = 


(4.2.1) 


Now,  given  (2.3.1)  together  with  the  boundary  condition  (2.3.2),  It  Is  easy  to  see 


that  the  following  holds: 
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where, 


^ 1  A  2,  A  3 

-'4  2  -^4  3  -4  4 


hv, 


^«rl 


An  A  i 


1 

01 

«1 

02 

u2 

0JVt 

UN1 

(4.2.2) 


2/i  =  (Xi, i.  ^.-,2.  •••.  y.-,w2)r 

and  «,•  =  (C/,-,!,  (/,-,2,  •  •  •  ,  Ui:NJT .  The  matrices  Au  A . .  are  all  N2xN2  clrculant 

matrices  depending  on  the  parameter  9  and  on  Q.  For  the  example  set  Q  of  sec¬ 
tion  (2.3) 


A i  —  o  -\/i  i,2 ,N1, 

A  j  =  Circulant  (1,03,0,.. .,0,^)  (4.2.3) 

A2=  Circulant  (04,06,O,...O) 

ANi  =  Circulant  (02,O,O,...O). 

The  linear  equation  system  (4.2.2)  has  a  block  circulant  matrix  structure  also,  l.e. 


where, 


B  ( 9)y  —  \fpu 


(4.2.2’) 


B  (9)  =  block  circulant  [A  X,A2,...,A^ 

Is  an  N1N2xNlN2  matrix. 

We  are  seeking  a  maximum  likelihood  estimate  for  9.  For  this  we  need  to 
write  an  expression  for  the  likelihood  function  for  the  field  {Y^}  or  equivalently 
the  random  vector  y.  Since  the  right  hand  side  of  (4.2.2’)  Is  an  NXN2  length 
Gaussian  vector,  It  follows  that  the  likelihood  function  Is, 


L  id, ft)  —  ±Lia  iB  n 


(ViftT™  ,s‘4i 

1  <j<n2 


exp{ - -(YiS  +  V  9kY.  2f} 

2/9  i  +  h  +h2'  ' 

Kashyap  [Kashyap  1980,  Image  Modeling]  has  shown  that, 


where, 


N,-l  n2-\ 

Det  (B  (0))  —  n  n  |  |  A  (4  ,4,0) 

1  =0  j  =0 
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(4.2.4) 


(4.2.5) 


*t1/-  A*2 


and 


A(2i>22,0)  *  l  +  Y,  M*ar(*i) 

A;  — 1 


2.-  =  «*P  [ - Tr - ]  *  =1.2. 


N, 


(4.2.6) 


(4.2.7) 


We  can  now  carry  out  a  maximum  (log)  likelihood  estimation  procedure  for 
the  unknown  parameters  (0,0)  using  the  finite  data  set  {Yitj  :  i <i  <N ,  i<j<N2}. 
Note, 


ln(L  (0,0)) 


w,  n2 


(N,N  2) 


ln(2  tt/9) 


E  E  lYa  +  E 


2d  ^ 

»  =1  J  =1 


k  =1 


(4.2.8) 


+  7  E  E  ln(|  \A(z\,zl,6)\  I2) 

^  I  =0  j  =0 


Maximizing  ln(L  (0,0))  with  respect  to  6  and  0  yields  the  estimates, 

0*  —  arg  [min  L  (0)] 

,  N2 

^  =-Njr,RBlY“  + 


k  =1 


+gk\j  +qk* 


(4.2.9) 


where, 
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Nl  N2 

L(0)  =  N1NM  E  E  (Yh. 

i=i  y=i 


+  E 


*=i 


Wrl  W2-l 

-  E  E  ln(  I  \A(z[,zi,0)  I  I2) 

1=0  y=o 

The  first  term  In  L  (9)  has  a  minimum  at, 

N !  N ,  JV2 

^  =  [  E  E  n?y)rr1  E  E  YtfYi, 

i  =i  y  =i  j  =i  y =i 

where, 


(4.2.10) 


(4.2.11) 


Yi>  "  +»? . ^  +^2)r  (4.2.12) 

The  formula  (4.2.11)  Is  very  useful  since  It  provides  a  quick  estimate  which  Is 

quite  close  to  the  optimal  6' .  In  fact,  one  may  use  a  Newton-Raphson  technique 

to  Iterate  from  1  to  6* . 

Everything  we  have  done  so  far  has  hinged  on  a  particular  choice  of  neigh¬ 
borhood  of  dependence  Q.  The  appropriate  choice  of  neighborhood  may  be 
decided  upon  by  considering  several  different  neighborhoods  Qv  Q2,...,  QR  and 
corresponding  to  each  Qk ,  compute 


Nr  1  N2- 1 

Ck  t  N XN 2  in  j3k  -  E  E  In  |  \  Ak(z[  ,z(  ,9k*)\  I  2  (4.2.13) 

*  ==o  y  =o 

where  0k,  9k  are  maximum  likelihood  estimates  of  the  pk ,  9k  In  the  model 
corresponding  to  Qk .  The  decision  rule  for  determining  the  neighborhood  of 
choice  Is 

select  that  Qk  which  minimizes  Gk ,  for  ke  {1,2 . R  }. 


See  [Chellappa  Kashyap  Ahuja  1979  Decision  rules  neighbors]  for  details. 
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4.3.  Segmentation/ Window  Determination: 

Our  main  goal  Is  to  describe  a  method  that  will  assist  us  In  Identifying  the 
window  size  NxxN2  over  which  one  can  expect  to  be  able  to  build  a  satisfactory 
homogeneous  model  for  the  gravity  and  bathymetry  fields.  Our  Idea  here  Is  that 
If  one  can  find  In  a  given  data  set  (gravity  or  bathymetry)  the  edges  across  which 
there  are  significant  changes  In  field  Intensity,  then  the  edges  can  be  used  as 
boundaries  of  segments/windows.  The  subject  of  segmentation/edge  detection  In 
Image  processing  has  a  very  large  literature.  See  the  book  [Rosenfeld  1981  Image 
Modeling]  for  diverse  points  of  view.  Kashyap  himself  suggests  a  statistical  tech¬ 
nique  based  on  hypothesis  testing  for  his  models.  We  find  that  his  methods  are 
rather  computation-intensive  and  we  present  an  alternative  originally  arising  In 
the  thesis  of  Canny  [Canny  1983  Edges  and  Lines].  The  method  Is  well- 
motivated  from  a  signal  processing  polnt-of-vlew  unlike  a  number  of  ad-hoc  pro¬ 
cedures  In  the  literature.  For  Instance,  It  brings  together  considerations  of 
slgnal-to-nolse  ratio  and  localization  Into  the  design  of  the  edge-detector.  It  leads 
to  algorithms  with  tunable  parameters  which  can  be  selected  to  adjust  detector 
performance.  It  Is  possible  to  motivate  the  main  Ideas  by  using  a  one  dimen¬ 
sional  model. 

4.3.1.  A  One-Dimensional  Edge  Model: 

Consider  a  one  dimensional  lmage/fleld  of  the  form, 

I  (x  )  —  Au  (x)  +  n(x) 

where  A  Is  a  constant,  n(x)  Is  noise  and  u(x)  is  the  unit  step: 


(4. 3.1.1) 
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( 1  x >0  , 

M(a;)={o  z  <0  (4.3. 1.2) 

Our  goal  Is  to  construct  a  procedure  for  detecting  the  edge  u  (•)  In  the  noisy  Image 
/(•)  and  to  localize  It.  We  consider  the  following  scheme: 


JlC*> 


#  £ 

Ooo 

JNJ  ■- 

*  i 

E  Cx) 


Figure  4.3.1 

Here  *f  denotes  convolution  by  a  suitable  kernel  function  f  and  N  denotes 
the  nonlinear  operation  of  suppressing  all  but  the  local  minimum  In  0(x).  The 
output  of  N  Is  supposed  to  be  the  edge/step.  The  first  (linear)  processor  Is 
chosen  to  be  a  convolution  by  requiring  translation  invariance  (l.e.  Insensitivity 
to  the  location  of  the  edge  -  here  the  point  x=0). 

Now, 

00 

0(0)  =  A  f  f(x)u(-x)dx  (4.3. 1.3) 

-  OO 

00 

+  f  f  (x  )n  (-x  )dx 

~  OO 

A  0  00 

=  jx  f  f  (x)dx  +  f  f  (x)n  (-x)dx  .  (4. 3. 1.4) 

-  00  ~co 

Assume  that  n(x)  Is  spatial  white  noise  with  zero  mean  and  “variance  parameter” 

oo 

<rn2.  Clearly,  0(0)  Is  a  random  variable  with  mean  —  A  f  f  (x)dx  and  the  spec- 

-  OO 


tral  energy  of  0(0)  (  or  O(x))  arising  from  the  noise  alone  Is  given  by, 
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E[  J  f  2(x  )n\-x  )dx  ] 

-  OO 

The  output  slgnal-to-nolse  ratio  SNR  Is  given  by 


(4.3. 1.5) 


A  f  f  (x  )dx 


S.N.R.  = 


f  f  2(x  )dx 

-  OO 

It  Is  desirable  to  choose  T  In  order  to  make,  SNR  as  large  as  possible. 


nr 


(4.3.1. 0) 


4.3.2  .  The  Localization  Problem: 


Since  the  nonmaximum  suppression  operation  marks  zero  crossings  of 


dx 


0  (x)  as  edge  points,  and  since, 


0=0  (x) 

=  o;  {x )  +  on'(x) 

OO 

=  A  }  {x)  +  f  f  (y  )n  (x  -y  )dy  ,  (4. 3. 2.1) 

-  OO 

It  can  be  shown  that,  under  the  assumption  that  f(x)  =  -f(-x),  the  quantity, 


A  I  f  (o) 


/'  \x)dx 


(4. 3. 2. 2) 


Is  a  measure  of  the  degree  of  localization.  If  L  Is  large  then  the  standard  devia¬ 
tion  of  the  distance  of  the  actual  maximum  to  the  true  edge  Is  small.  Consider, 


J  /  (x  )dx 


S  =  £(/  )  = 


A7 


/  z(x  )dx 


(4. 3. 2. 3) 


Both  E  and  A  are  Independent  of  the  magnitude  scales  for  f.  Furthermore, 
suppose  we  do  spatial  scaling, 


/()  -  /*(•) 

where, 


fw(x)=f  (x/w). 

Then,  It  can  be  verified  that, 


=  W  )  Mf  ) 


Thus  localization  can  be  traded  off  against  detection  performance.  A  ‘broad’ 
function  /  (•)  will  have  good  signal  to  noise  ratio  but  poor  localization  compared 
to  a  filter  with  a  narrow  /  (•).  One  could  now  try  to  find  /  (•)  that  maximizes 
E (/  )A (/  ).  One  can  also,  for  computational  ease,  try  to  find  f  that  maximizes 
A (/  )  while  E (/  )  =  Ci  Is  fixed  or  find  f  that  maximizes  E (/  )  while  A (/  )  =  c2  Is 
fixed.  The  corresponding  Euler-  Lagrange  equations  yield  explicit  families  of  /  (•) 
that  are  appropriate  as  filters  [Canny  1983  Edges  Lines]. 
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The  principal  defect  of  the  above  approach  Is  that  It  does  not  take  Into 
account  the  possibility  multiple  edge  responses  from  the  detector  caused  by  noise. 
From  the  basic  edge  detection  scheme  (see  figure  4.3.1),  It  Is  clear  that  such  mul¬ 
tiple  detector  responses  are  due  to  zero  crossings  In  the  derivative  of  the  output 
of  the  linear  filter.  These  (random)  zero  crossings  can  be  estimated  by  a  formula 
due  to  S.  Rice  and  we  are  lead  to  an  optimization  problem  of  the  form, 

0 

minimize  /  (/  2  +  X,/'  2  +  X2/'  '  2  +  \J  )dx  (4.3.2. 6) 

-w 

where  Xj,  X2,  X3,  are  Lagrange  Multipiers. 

Analysis  of  the  corresponding  Euler-La^grange  equations  leads  to  a  family 
of  solutions, 

/  {x)  —  die"  sin(ojx)  +  a  2e  “  cos  (ojx  ) 

+  a  3e~ax  sin(wa: ) 

I 

+  a  ie~ax  cos  (ujx)  +  c 

with  boundary  conditions, 

/  (0)  =  o  =  /  (-W) 
f  (0 )  =  s;  /'  (-Vk)  =  0. 

These  lead  to  conditions  on  the  parameters  a,-  and  c.  In  terms  of  s.  We  omit  the 
details. 

The  nonlinear  processor  N  for  nonmaximum  suppression  Is  quite  standard  In 
the  literature  on  edge  detection/segmentation  and  we  refer  to  page  81  of  [Canny 
1983  Edges  and  Lines]. 
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We  have  Implemented  a  discrete  2-D  version  of  the  above  edge  detection 
algorithm  Involving  the  following  steps. 

Let  /(»,/)  be  the  given  Image.  Let  {A(/)}  the  discrete  convolutor. 

Step  1:  Compute  Y(1,J), 

oo 

Yr(i.j)=  E  A  (/-/)/(*,/) 

l  —1  -  OO 

Step  2:  Compute  X(1J) 

OO 

X(i.j)*=  E  A  (i-l)YU.j) 

l  =  -  OO 

Step  3:  Do  nonmaximum  suppression.  The  above  Infinite  convolutions  are 
truncated;  since  we  use, 

hk  =  e'x*  cos  (wk  )  k  =0,1,2,... 

A_*  = -A*  A  =0,1,2,... 

It  Is  possible  to  Implement  the  above  algorithm  using  /as£  convolutions.  The 
parameter, 

.  l 

X  =  —  >  o, 

er 

and 

w  =  0.8X. 

Tuning  cr  gives  us  control  of  performance.  A  smaller  value  of  <r  give  us  good 
localization,  while  a  bigger  value  of  o  picks  up  global  features.  A  series  of  com¬ 
puter  experiments  were  performed  on  a  VAX  11/780  computer  at  the  Computer 
Vision  Laboratory  at  the  University  of  Maryland,  using  digitized  pictures  of  real 


scenes  and  objects.  By  tuning  the  a  parameter  appropriately  It  was  possible  to 


(c) 


Figure  4.4.1 

(a)  object 

(b)  Canny  detector 

*  =  -V3.45 
co  =  0.7 

(c)  Sobel  mask 


(e) 


Figure  4.4.1 

(d)  urban  scene 

(e)  Canny  detector,  parameters 
same  as  in  figure  4.4.1  (b) 

(f)  Sobel  mask 


Figure  4.4.1 

(g)  portrait 

(h)  Canny  detector 

*  =  +1.0 

Co  =  0.7 


Note:  all  pictures  with  255  gray  levels, 
image  sizes  128  x  128  pixels. 
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often  obtain  better  performance  In  edge  detection  than  with  the  help  of  a  tradi¬ 
tional  algorithm  such  as  the  one  based  on  a  differencing  mask  due  to  Sobel. 
Some  typical  pictures  are  Included  (see  Figure  4.4.1). 

4.4.  Some  Recent  Advances  in  Deconvolution: 

Here  we  note  rather  briefly  how  some  new  developments  In  deconvolution 
may  be  useful  for  geophysical  Inversion. 

Many  observation  models  In  geophysical  exploration  take  the  form, 

/  *  A*  -  /  /  (*  -  y)dy(y)  =  g,  (x)  t=i,2,..JV  .  (4  4i) 

R2  ‘  V  ‘  ' 

where  p  Is  a  distribution  of  compact  support  over  a  region  In  R”2  defining  an 

t 

observation  process  and  /  (•)  denotes  an  unknown  geophysical  field  to  be  deter¬ 
mined  (e.g.  density  anomaly)  and  &•(•)  denotes  the  corresponding  observation 
(gravity  anomaly,  magnetic  anomaly  etc.).  The  Indexing  of  the  observation  data 
set  by  i  e{i,2,...,N  }  may  be  used  also  to  refer  to  multiple  passes/tracks. 

The  Inverse  problem  of  Interest  to  us  Is  to  determine  /  (•)  given 
£,•(•)>  i  E{i,2,...,N  }.  We  call  this  a  deconvolution  problem  involving  multisensors.  If 
however,  only  one  data  set  say  g1  Is  available,  then  It  Is  well-known  that  the 
Inverse  proble  Is  Ill-posed  --  the  Inverse  opertor  Is  typically  unbounded  —  and 
leads  to  serious  Issues  of  numerical  accuracy.  Various  regularization  techniques 
are  known  for  the  problem  of  approximate  reconstruction  of  /  (•). 

On  the  other  hand,  If  we  have  multiple  sensors,  and  If  we  can  solve  for  o,-  In 


the  equation, 
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ul  *  fj,  +  *  fi  +  +  i /N  *  n  =  8  (4.4.2) 

where  S  denotes  the  Dirac  delta  function,  then  It  Is  clear  that  we  have  a  recon¬ 

struction  formula, 

/  =  *  9i  +  v-i  *  9i  +  '  '  ’  +  uN  *  gN  •  (4.4.3) 

By  taking  Fourier  transforms  of  both  sides  In  (4.4.2),  we  obtain  the 

equivalent  form, 

W  +  v2/i  +  ■■■  +  vNn  =1  .  (4.4.4) 

Equation  (4.4.4)  Is  known  as  the  Bezout  equation.  In  their  fundamental 
work,  Berensteln,  Taylor  and  Yger  obtained  explicit  formulas  for  the  deconvolu- 
tors  Vi,...,vN,  for  certain  classes  of  problems.  The  key  requirement  here  Is  that,  In 
order  for  the  existence  of  the  deconvolutors  vu  .  .  .  ,  uN  (equivalently,  solvability 
of  the  Bezout  equation),  the  functions  n  ,...n  should  satisfy  a  coprlmeness  condl- 

l  N 

tlon,  l.e.,  In  a  strong  sense,  these  Fourier  transforms  should  not  have  common 
zeros. 

In  a  long  report  [Berensteln  Krlshnaprasad  Taylor  1985  Deconvolution 
methods],  the  authors  undertook  a  detailed  numerical  study  of  a  special  class  of 
examples.  These  examples  Involved  two  sensors  (N=2)  and  the  observation  map 
corresponded  to  averaging  over  an  Interval  »=i,2,;  thus 

x  +  at 

9i(.x)  =  ~-  /  f  (y)dy  i=i,2  . 

2a>  *-a, 

the  corresponding  Fourier  transforms  n  and  n  are  given  by 

1  2 


(4.4.5) 
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sin(  a,-  x ) 

V  =  — — -  *'  =1,2 


(4.4.6) 


The  necessary  and  sufficient  condition  for  strong  coprlmeness  (solvability  of 
the  Bezout  equation)  turns  out  to  be: 


ax/a2  is  irrational  !!  (4.4.7) 

In  the  paper  [Berensteln  Krlshnaprasad  Taylor  1985]  extensive  numerical 

studies  were  conducted,  using  aj=i,  a2=\/2.  The  technique^  applied  to  the  prob¬ 
lem  of  resolving  peaks  In  a  double  Gaussian  proved  to  be  remarkably  robust. 

Our  main  point  here  Is  that  It  would  be  very  useful  to  try  and  attempt  a 
study  of  the  suitability  of  the  above  deconvolution  methods  to  problems  of  geo¬ 
physical  Inversion.  As  a  first  step  we  suggest  the  use  of  this  method  on  transfer 
function  models  relating  bathymetry  and  free  air  gravity. 


5.  Future  Work: 

In  this  section  we  make  some  preliminary  recommendations  for  further 
development  of  our  efforts  In  order  to  provide  practical  tools  for  basic  geophysical 
modeling  and  Inversion  techniques.  The  pervasive  theme  of  this  report  Is  that 
while  many  successful  algorithms  have  been  produced  for  these  problems,  these 
tend  to  be  primarily  based  on  1-D  processing.  With  the  recent  advances  In  the 
subject  of  discrete  random  fields,  It  Is  desirable  to  take  advantage  of  new  2-D  sig¬ 
nal  processing  methods.  There  are  several  possible  directions.  As  a  first  step.  It 
Is  desirable  to  test  the  proposed  2-D  statistical  modeling  and  Inversion  techniques 
on  synthetic  and  real  data.  Detailed  sea-bottom  models  should  be  Investigated 


with  reference  to  their  correlation  properties. 
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5.1.  Testing  of  Proposed  Inversion  Technique: 

We  suggest  that  the  modeling  and  Inversion  techniques  of  section  4  be  Imple¬ 
mented  and  first  tested  on  synthetic  data  of  reasonable  verisimilitude.  Special 
attention  should  be  given  to  the  types  of  parameter  ranges  and  Q  neighborhoods 
Involved.  The  segmentation  algorithm  should  also  be  Implemented  and  tested 
alongside. 

It  would  be  very  natural  to  treat  the  2-vector  ( G{]- ,  Bi}  )T  as  a  bivariate  ran¬ 
dom  field  and  then  proceed  to  construct  models  for  It  directly  using  data  from 
regions  where  we  know  both  gravity  and  bathymetry.  However,  at  present  this 
does  not  appear  to  be  computationally  feasible.  It  Is  for  this  reason  that  we  have 
chosen  to  work  with  univariate  field  models.  However,  It  Is  desirable  that  such 
bivariate  modeling  methods  be  further  developed  and  fast  algorithms  should  be 
Investigated. 

An  Important  part  of  our  methodology  Involves  the  modeling  of  bathymetry. 
It  Is  desirable  that  detailed  statistical  studies  of  sea-bottom  roughness  be  con¬ 
ducted  In  order  to  get  some  insight  Into  this  basic  element  of  our  methodology. 
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